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We propose a method for Monte Carlo simulation of sta- 
tistical physical models with discretized energy. The method 
is based on several ideas including the cluster algorithm, the 
multicanonical Monte Carlo method and its acceleration pro- 
posed recently by Wang and Landau. As in the multibondic 
ensemble method proposed by Janke and Kappler, the present 
algorithm performs a random walk in the space of the bond 
population to yield the state density as a function of the bond 
number. A test on the Ising model shows that the number 
of Monte Carlo sweeps required of the present method for 
obtaining the density of state with a given accuracy is pro- 
portional to the system size, whereas it is proportional to the 
system size squared for other conventional methods. In ad- 
dition, the new method shows a better performance than the 
original Wang-Landau method in measurement of physical 
quantities. 



I. INTRODUCTION 

The Monte Carlo simulation is one of the most pow- 
erful tools for investigating models in statistical physics 
Although the Metropolis method Q and its varia- 
tions are available for simulating variety of models, they 
are not necessarily the best methods when the system of 
interest has a strong long-ranged correlation. Essentially 
two approaches have been proposed for overcoming the 
drawbacks of such local-updating methods. In one ap- 
proach, one uses an ensemble entirely different from the 
ordinary canonical ensemble with a fixed temperature, 
whereas in the other approach one extends the original 
ensemble by introducing auxiliary variables. 

Multicanonical method (|-[|] , broad histogram method 
H , and the flat histogram method JtJ belong to the first 
category. In these methods a random walk in the en- 
ergy space is performed to calculate the state density 
as a function of the energy. Multicanonical method was 
applied to the Q-state Potts model, for example, and 
turned out very successful [||. Meanwhile it was real- 
ized that the random walker tends to be blocked by the 
edge of the already visited area. In addition, because of 
the general feature of random walks, it takes a long time 
to go from one end of the area to the other. Recently 
Wang and Landau Q succeeded in removing these prob- 
lems by penalizing moving to and staying at the energy 
which has been visited many times. The efficiency of the 
Wang-Landau (WL) method was also demonstrated in 
an application to antiferromagnetic Q-state Potts model 
on a simple cubic lattice ||. In particular, the method 
turned out to be powerful in studying the ground state 



properties due to the fast diffusion accelerated by the WL 
method. 

The second category includes various cluster algo- 
rithms. In cluster algorithms, graph degrees of freedom 
are introduced to extend the original ensemble. In most 
of their successful applications, clusters of the size of cor- 
relation length are formed and flipped. A cluster 
algorithm is applied to Q-state Potts model and proven 
much more efficient f!o|| than local updating algorithms. 

Janke and Kappler |T^] proposed the multibondic 
method, the combination of the multicanonical method 
and the cluster algorithm. They took the two- 
dimensional Q-state Potts model as an example. They 
measured the computational time required for the ran- 
dom walker to traverse the interval between the two 
peaks in the canonical probability distribution. The tra- 
verse time in units of Monte Carlo sweeps was found to 
be proportional to the number of spins, N , whereas it 
is proportional to N for the ordinary multicanonical 
method, where N is the total number of spins in the 
system. The multicanonical cluster algorithm Jl3| shows 
the same size dependence as the Janke-Kappler algorithm 
(referred to as 'JK' in the present paper). The compari- 
son between the Janke-Kappler algorithm and this algo- 
rithm can be found in p"2| . 

In this paper, we propose a new method based on the 
JK method (or, more generally, the multibondic ensem- 
ble method) and the Wang-Landau acceleration method. 
In Section II, we briefly review the multibondic ensem- 
ble method. In Section III, we propose a modification of 
the Janke-Kappler method (MJK) to avoid the possible 
lasting effect of the initial graph. Then, in Section IV, a 
combination of the MJK with the Wang-Landau method 
(WL) is discussed. We refer to this combined method 
as MJKWL. In Section V, we demonstrate the efficiency 
of the MJKWL by comparing it with other methods. In 
particular, it is shown that the MJKWL is better than 
each of its ingredients, i.e., the JK method and the WL 
method. In Appendix A, we present a simple and exact 
relationship between the density of state (DOS) as a func- 
tion of the energy and the DOS as a function of the bond 
number for the Q-state Potts model in any dimensions. 

Since we are forced to use many abbreviations in the 
present paper to refer to various methods, it may be con- 
venient to summarize all of them here: 

JK ... The Janke-Kappler method. 

WL ... The Wang-Landau method. 

SWL ... The Wang-Landau method with single spin up- 
date, i.e., the original Wang-Landau method. 
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MJK ... The modified Janke-Kappler method. 

MJKWL ... The modified Janke-Kappler method with 
the Wang-Landau acceleration method. 

II. THE MULTIBONDIC ENSEMBLE METHOD 

Since our method can be viewed as a derivative of the 
Janke-Kappler algorithm, we give a brief review of the 
multibondic ensemble method. In the following, we take 
the Q-state Potts model as an example to make our de- 
scription concrete. However, the generalization to other 
models with discrete energy is straightforward, and we 
try to describe our method so that the generalization ap- 
pears obvious. 

The Hamiltonian is given by 

H = -^X, <7i,CT 3 > °i = {L • • • lQ} 

(ij) 

where J is the exchange coupling constants and (ij) de- 
notes a nearest neighbor pair. In what follows, we take J 
as the unit of the energy, and J / Kp as the unit of tem- 
perature, where kp is Boltzmann's constant. We first 
represent the partition function as a double summation 
over states S and graphs G, following the general frame- 
work of the dual algorithm ; 

Z(T) = W (S) = W (S, G)=J2 MG) A(S, G) 

S S,G S,G 

(1) 

This is nothing but the well-known Fortuin-Kasteleyn 
representation ]l7|]. Wo (S) is the weight of state 5, 
whereas A(S, G) is a function that takes the value one 
when S is compatible to G and takes the value zero oth- 
erwise. Vq(G) denotes the weight for graph G defined 
as 

Vo(G) = V (n b (G),T) = (e l l T - 1)"^ 

where n b (G) is the number of bonds in G, in the case 
of the Q-state Potts model. Although the only graph el- 
ements are bonds in this case, a graph consists of more 
than one type of elements in general applications. There- 
fore, in more general terms, n b (G) is a p-dimensional vec- 
tor variable i-th element is the number of graph elements 
of the i-th kind contained in the graph G. By taking the 
summation over S and G, fixing the fixed number of the 
bond, the above expression for the partition function is 
reduced to 

Np 

Z(T) = Q(n b )V (n b ,T). (2) 

n b =0 

where Np is the total number of nearest neighbor pairs 
in the whole system (Np = dN = dL d for d-dimensional 



hyper cubic lattices). Here, fl(n b ) is the DOS of the bond 
number defined as the number of consistent combinations 
of graphs and states such that the graph consists of n b 
bonds; 

n(n b )= Yl E A ^ G )- 

{G\n b (G)=n b } S 

In the multibondic ensemble method, we replace 
Vo(n b ,T) by a free function, which we denote by V(n b ). 
By adjusting this function, we try to make the histogram 
flat as a function of n b . In other words, V(n b ) is adjusted 
so that the product of fl(n b ) and V(n b ) may be indepen- 
dent of n b . Since initially Q(n b ) is unknown, we achieve 
this adaptively. The way how we modify the trial weight 
is described in the following sections. 

There are several different ways for adjusting the trial 
weight. The original JK method is one of them. The 
detailed description of the JK method can be found in 
Appendix B and in the original paper jl2| |. 

III. MODIFICATION OF THE JANKE-KAPPLER 
ALGORITHM 

As stated in Appendix B, in the original JK algorithm 
the initial graph for each Monte Carlo sweep may have a 
long lasting effect upon the subsequent states during the 
sweep since the bond update is done only sequentially. 
In this section, we propose a modification of the Janke- 
Kappler algorithm (MJK) to reduce the possible lasting 
effect of the initial graph as much as possible. 

We first choose any consistent combination of a state 
and a graph as an initial condition. For the initial choice 
of V(n b ), we choose V(n b ) = 1 for any n b . 

In each Monte Carlo sweep of the MJK, we choose the 
number of the bonds to be placed on the whole system 
before actually placing them. The number n b is chosen 
with the following probability, 

Here, n p (S) is the number of "satisfied" pairs in the cur- 
rent state S, i.e., 

n p( S ) = J2^ {S ) :<J] (.S)- 

iv) 

and 

(i\ s _a_ 

\ m J m \(l - m )l 

Note that the choice of n b is based only on the informa- 
tion of the current state S at the beginning of the sweep. 
As for the graph G at the beginning of the sweep, we sim- 
ply delete all the bonds in it. We then choose n b pairs at 
random out of n p (S) satisfied ones and place new bonds 



2 



on them. It is clear that there is no direct influence of the 
initial graph on the final graph. The correlation between 
them arises only through the state which should be com- 
patible to the initial graph. After the placement of the 
n b bonds, we 'flip' all the clusters of sites with probabil- 
ity one half where flipping a cluster means changing all 
the variables on it simultaneously. This completes one 
Monte Carlo 'sweep'. It is easy to show that this pro- 
cedure satisfies the extended detailed balance condition 
& 

P{G\S)W{S) = P(S\G)W(G) 

where W(G) = T,s W ( S > G ) and W(S,G) = 
V{G)A(S, G). 

Another modification should be done in order to make 
MJK method better than the original JK method. In 
the JK method, the histogram H{n b ) is updated by the 
simple rule 

H{n b )^H(n b )+S(n b ,n b (G)) (3) 

every time a part of graph, i.e., a bond on a pair of sites, 
is updated. It means that we update H(n b ) in a sweep as 
many times as the number of bonds. Although the suc- 
cessive values of n b in the same sweep are strongly corre- 
lated with each other in the JK algorithm, one can still 
get statistically more informative data by taking them 
all into account. It is roughly equivalent to adding some 
smooth function to H{n b ) at every Monte Carlo sweep, 
in contrast to adding a delta function. 

However, in the MJK method, H{n b ) is updated only 
once in every Monte Carlo sweep, which means that a 
delta function is added to H(n b ) at each sweep according 
to the the updating rule Eq. (|J) . In order to remove this 
disadvantage, in the MJKWL method, we add to the 
histogram the expectation values of the delta function 
5(n b ,n g (S)), rather than the delta function itself. The 
resulting updating rule for H(n b ) is 

H{n b ) <= H{n b ) + N P P(n b \n p (S)). 

Although including a constant Np is not relevant, it is 
added in order to make the magnitude of the histogram 
comparable to the one in the JK method. 

IV. COMBINING WITH THE WANG-LANDAU 
METHOD 

It has been known since the first proposal of the mul- 
ticanonical method that the random walker tends to be 
stuck at the boundary which separates the region visited 
already from the one not visited yet. This difficulty has 
been removed by the recent technique proposed by Wang 
and Landau Q . Their method seems to be useful also in 
accelerating the diffusion of the random walker. 

In the WL method, the DOS phase (see the next sec- 
tion) of the computation consists of varying number of 



consecutive sets of simulation, as is also the case with 
the ordinary multi-canonical method and its derivatives. 
However, the important difference lies in the way the trial 
weight is updated. In the conventional multi-canonical 
methods, the trial weight is updated only at the end of 
each set of simulations. During each set, the trial weight 
and, consequently, the transition probability are fixed. 
In the WL method, on the other hand, every time the 
state of the system is renewed, the trial weight W(E) is 
updated as 

hiW(E) <J= InW(E) - X5(E,E(S)) (A > 0), 

where E(S) is the energy of the current state. The posi- 
tive parameter A is introduced to control the magnitude 
of the expelling force imposed on the random walker. If 
the parameter is large, the random walker quickly moves 
out of the region which it has already visited. In other 
words, the histogram is forced to be flat by this param- 
eter. However the very presence of this force breaks the 
detailed balance and therefore makes the resulting trial 
weight differ from the correct one, i.e., the inverse of 
the DOS. On the other hand, if it is small, the resulting 
trial weight is reliable while the convergence tends to be 
slow. Therefore, the general strategy is to set this value 
large initially and make it smaller as the trial weight ap- 
proaches the correct one. 

The value A is kept constant throughout each set of 
simulation and is reduced by some factor at the begin- 
ning of the next set of simulation. Wang and Landau 
suggested 1/2 for the reduction factor. Each set of sim- 
ulation terminates when the histogram of the set sat- 
isfies some predetermined condition concerning its flat- 
ness. The histogram is reset at the beginning of a new 
set while the trial weight W(E) is not. The whole cal- 
culation is terminated when A becomes smaller than a 
predetermined value. 

In order to combine the WL method with the MJK 
method described in the previous section, we may simply 
replace the energy space in the original WL method by 
the bond-number space. To be more specific, H(E) is 
replaced by H(n b ), W{E) by V(n b ), and S(E,E(S)) by 
6(n b , n b {G)). However, as for the updating rule of V(n b ), 

]nV(n b ) <= bxV(n b ) - XN P P(n b \n p (S)) 

is the better choice than 

\nV{n b ) <= \nV(n b ) - \6{n b ,n b (G)) 

for the same reason as we stated in the previous section. 
We therefore use the former updating rule in the MJKWL 
method for the sample calculation presented below. 
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V. EFFICIENCY OF THE METHOD 

We now discuss the performance of the above- 
mentioned methods. Since the present method 
(MJKWL) is the combination of MJK and WL, it should 
be demonstrated that this combination is meaningful, 
i.e., MJKWL is qualitatively better than both of the two 
ingredients. 

First it should be noted that there are several measures 
of performance. In the WL method, such as MJKWL 
and SWL, the whole computation process consists of two 
phases; the DOS phase and the measurement phase. In 
the DOS phase, the computation is performed mainly 
for obtaining an estimate of the DOS. During this phase, 
several sets of simulation are done for the adoptive ad- 
justment of the fictitious weight (the trial weight V(tib) 
in Section || for the multibondic methods such as JK, 
MJK and MJKWL and the trial weight W(E) in Section 
[V for the Wang-Landau method with single spin update 
(SWL)). At the end of this phase, some of the physi- 
cal quantities, such as the entropy, the energy and the 
specific heat, can be computed with the resulting DOS. 
For other quantities, however, some additional simula- 
tion should be performed with a fixed trial weight and 
with the controlling parameter A set to be zero. We call 
this part the measurement phase. (In the case where 
the value of A for the last set of simulation in the DOS 
phase is negligibly small, a separate measurement phase 
may not be necessary. In such cases, we regard the last 
set as the measurement phase.) In what follows, we dis- 
cuss the computational time required for the DOS phase 
and that for the measurement phase, separately. As we 
see below, the DOS is obtained much faster in the DOS 
phase of MJKWL than in JK and MJK, while SWL shows 
qualitatively the same performance as MJKWL. The dif- 
ference between MJKWL and SWL can be seen in the 
measurement phase. Namely, MJKWL yields much bet- 
ter statistics in the measurement phase than SWL within 
the same number of Monte Carlo sweeps. 

A remark should be placed here concerning the sources 
of errors in the two phases. During the DOS phase, the 
systematic error as well as the statistical error is present. 
The systematic error is due to the obvious fact that there 
may be a region which the random walker has not visited 
yet. In the methods based on the Wang-Landau accel- 
eration, the fact that the controlling parameter A is not 
zero in another source of systematic error. Because of 
this systematic error, the dependence of the total error 
on the duration of the simulation is complicated. On 
the other hand, in the measurement phase there is no 
other sources of errors than the ordinary statistical ones. 
Therefore the precision of the result in this phase is pro- 
portional to the inverse of the square root of the number 
of Monte Carlo sweeps. 

In what follows, we argue and demonstrate that the 
methods without the Wang-Landau acceleration, such as 
JK and MJK, require the number of Monte Carlo sweeps 



of the order 0(N 2 ) whereas the methods with the Wang- 
Landau acceleration, such as MJKWL and SWL require 
O(N) to achieve the same accuracy in the DOS. 

In all the methods, discussed here, we start with some 
ad-hoc initial guess for the DOS. Then, the resulting his- 
togram has a rather narrow range of distribution. There- 
fore, in order to make the histogram flat throughout the 
whole energy (or bond number) range, we have to repeat 
simulations. Every time we start a new set of simulation, 
we improve the initial guess for the density of states based 
upon the outcome of the last set of simulation. The diffi- 
culty arises near the boundary between the two regions; 
the region which has been visited already in the previous 
simulations and the region which has not. When the ran- 
dom walker in the energy (or the bond number) space hits 
the boundary during the simulation, it usually bounces 
back and, even if it does not, it seldom goes far beyond 
the boundary. Therefore, the width of the visited region 
increases by only a few steps as a result of the whole set. 
It follows that the number of sets of simulation required 
for making the histogram flat is proportional to the width 
of the energy (or bond number) space, that is, O(N). In 
addition, each set must be long enough for the walker 
to traverse the whole previously visited region. Since the 
width of the previously visited region is of the order O(N) 
in general and the typical distance the walker traverse in 
a single Monte Carlo sweep is OiN 1 / 2 ), the number of 
Monte Carlo sweeps required for the walker to traverse 
the region is 0({N /N 1 / 2 ) 2 ) = O(N). These factors are 
multiplied to make the total number of sweeps required 
for the whole DOS phase of the order 0(N 2 ). 

In contrast, in methods with Wang and Landau's ac- 
celeration, the situation described above cannot hap- 
pen. This is because the current histogram affects the 
current weights and transition probabilities, such that 
the weights for the frequently visited positions become 
smaller. This forces the walker to move out of the al- 
ready visited region and make histogram flat. Since this 
situation is similar to the one-dimensional self-avoiding 
walk, a natural guess is that the number of steps (i.e., 
the number of local updatings) required for the walker 
to traverse the already visited region is proportional to 
the size of the region, which is O(N). In units of sweeps, 
it is 0(N°), 

To check if this simple argument is correct, we per- 
formed simulations for ferromagnetic Ising model on a 
square lattice using three different methods: MJKWL, 
MJK and JK. For these three methods, we set an initial 
weight V( n f>) = 1 f° r au n b- We measured the number 
of Monte Carlo sweeps required for obtaining the DOS 
with a roughly fixed precision, as a function of system 
size N = L 2 . It should be remarked here that we cannot 
rigidly fix the target precision of the DOS because the 
termination condition in MJKWL is defined in terms of 
the flatness of the histogram and the value of the control- 
ling parameter A, not the number of Monte Carlo sweeps 
nor the precision of the DOS. Therefore, we performed 
a MJKWL simulation first with some reasonable choice 
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of the termination condition. Then, we performed simu- 
lations using JK and MJK. The current estimate of the 
DOS is updated frequently in these simulations so that 
the simulation can be terminated as soon as the precision 
of the DOS estimate reaches the same as that obtained 
in the MJKWL simulation. The precision of the DOS is 
measured by the following quantity. 
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- M(»6) -ln^ cxact >(n b ) 



The exact DOS, ft( cxact ), is obtained through Eq. (§) 
and the exact DOS as a function of the energy [|4|. In 
what follows, the termination condition for the MJKWL 
is the same for all the system sizes. It turned out that the 
resulting signal-noise ratio of the DOS, e(L), is roughly 
independent of the system size. 

Our procedure for the MJKWL simulation is as follows. 
The reduction factor A is divided by two when each set 
of simulation is terminated. Each set is terminated when 
the smallest H(rib) becomes greater than 0.8 times the 
average value of H(rib). The whole calculation is termi- 
nated when A becomes less than 10~ 8 . This procedure 
is essentially the same as suggested in the original paper 
by Wang and Landau (?[ except that we work with the 
DOS as a function of nj, rather than E. 

For MJK and JK, we perform a number of subsequent 
sets of simulations to improve the estimates of the DOS. 
We start with a relatively short set and gradually make 
it longer. The way we increase the number of sweeps 
of a set depends upon whether the random walker has 
already visited the whole bond-number space. If rib has 
not visited the whole rib space at the end of the 7-th set, 
the number of sweeps for the i + 1-th set is chosen as 

U+i = 10(771* + y/Np) 

where m, is the number of the already visited values of 
jib in the 7-th set. If the above argument is correct, i.e., 
the random walker moves in the bond-number space as 
a self-avoiding walker, this choice of ti+\ should give the 
walker an enough time to traverse the whole region of 
previously visited values of rib and touch the boundary 
a few times. Therefore it should be enough to expand 
the visited region. If rib has already visited the whole rib 
space at the end of the 7-th set, the number of sweeps for 
the i + 1-th set is given by 



t 



i+i 



= 2U. 



This choice will provide the walker with an enough time 
to develop appreciably better trial weights than the pre- 
vious sets. The whole procedure yields the total number 
of sweeps of the order 0{N 2 ) if the above argument is 
correct. The entire process is terminated when the esti- 
mated n(nb) has become as good as that obtained with 
the MJKWL. 



The computation is done for system sizes L = 
4, 8, 16, 24 and 32 as all other sample calculations pre- 
sented below. The results is shown in Figure 1. The 
average is taken over about 30 independent simulations. 
We can easily see that MJKWL is the best method among 
the three multi-bondic methods for larger N. It can be 
also seen that the MJK is better than the JK. Two lines 
are drawn in Figure 1 for references. The lower dashed 
line corresponds to t oc 0(N 1 ) whereas the upper dashed 
line t oc 0(N 2 ) We can see that the MJKWL requires 
0(N 1 ) sweeps while the MJK and the JK require 0(N 2 ) 
sweeps, as expected from the argument. We have also 
confirmed that the relative statistical error in the DOS 
obtained by MJKWL does not strongly depend on the 
system size. 
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FIG. 1. The total number of Monte Carlo sweeps per- 
formed to obtain the same accuracy in the DOS estimate, 
as a function of the system size N = L 2 for ferromagnetic 
Ising model on a square lattice. Three methods are examined: 
the modified Janke-Kappler algorithm with the Wang-Landau 
method (MJKWL), the modified Janke-Kappler algorithm 
(MJK), and the ordinary Janke-Kappler algorithm (JK). The 
average is taken over 30 independent runs. The upper and 
lower lines are for references, corresponding to t oc 0(N 2 ) 
and t oc 0(N 1 ), respectively. 



The performance of SWL, i.e., Wang and Landau's 
original method using single spin-flips, is also examined. 
We set the initial weight W(E) = 1 for all E. We mea- 
sured the total number of Monte Carlo sweeps as a func- 
tion of the system size. Again the precision of the result- 
ing estimate of the DOS does not strongly depend on the 
system size. The result is shown in Figure 2. We can 
see that the SWL require 0(N 1 ) sweeps for large sys- 
tems. Therefore, it can be concluded that the SWL has 
the same qualitative performance as the MJKWL in the 
DOS phase. 
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FIG. 2. The total number of Monte Carlo sweeps as a 
function of the system size N = L 2 for ferromagnetic Ising 
model on a square lattice. The calculation is performed fol- 
lowing the Wang and Landau's original procedure. The av- 
erage is taken over 10 independent runs. The upper and 
lower lines are for references, corresponding to f a 0(N 2 ) 
and t oc 0(N ), respectively. 
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FIG. 3. The standard deviation of 50 independent esti- 
mates of the squared magnetization per spin thermally aver- 
aged at the critical temperature, Tc = 2/ln(l + y2), for the 
ferromagnetic Ising model on a square lattice. For each run, 
100 x N sweeps are performed. 



To compare the efficiency of MJKWL and SWL in 
the measurement phase, we calculate the magnetization 
squared, M 2 divided by N 2 for ferromagnetic Ising model 
on a square lattice. We first estimate the DOS in the 
DOS phase. Using this DOS, we then perform 50 inde- 
pendent runs for the measurement phase using different 
random number sequence for each run. Each run con- 
sists of 100 x N sweeps, and produces a histogram and 
a set of microcanonical averages of the squared magne- 
tization, as is usually done in any multicanonical-type 
methods. Based on these information, the canonical av- 
erage of the squared magnetization at the critical temper- 
ature is computed for each run. Then, we compute the 
standard deviation of these 50 canonical averages. This 
standard deviation is proportional to the statistical error 
in the final estimate and can be used as a measure of the 
efficiency with which the spin configuration is updated 
during the simulation. 

The result is shown in Figure 3. As is clear from the 
figure, MJKWL is better than SWL. The difference in 
the standard deviation tends to increase as the system 
becomes larger. This is because the spin configuration 
is updated by clusters in MJKWL whereas it is updated 
by single spins in SWL. Therefore, the configuration is 
decorrelated much faster in MJKWL than in SWL. To 
be more specific, a random walker in SWL must visit 
states with very different values of energy in order to visit 
a state with very different value of the magnetization, 
whereas a random walker in MJKWL does not have to 
because the state can change even without changing the 
bond number at all. 



VI. SUMMARY 

We have proposed a combination of the Janke-Kappler 
algorithm with the Wang-Landau acceleration method, 
together with a modification of the Janke-Kappler algo- 
rithm. The number of Monte Carlo sweeps required for 
obtaining the DOS with several methods have been mea- 
sured and compared. It has been demonstrated that the 
number of Monte Carlo sweeps required for obtaining the 
DOS in the methods without the Wang-Landau accelera- 
tion is proportional to iV 2 , whereas in the present method 
(MJKWL) it is proportional to TV 1 . The new method is 
also compared with Wang and Landau's original method 
based on single spin updating (SWL). The result shows 
that the spin configuration is much more efficiently up- 
dated in MJKWL than in the SWL. 

The proposed modification to the Janke-Kappler algo- 
rithm turns out to be useful in reducing the CPU time 
requirement, though not as vital as the Wang and Lan- 
dau's idea in the cases shown in the present article. 

We have also deduced an exact relation between the 
DOS as a function of the energy and that as a function 
of the bond number for the Q-state Potts model in any 
dimensions. 

The present method can be easily extended to other 
models with discrete degrees of freedom, in particular, 
when the cluster algorithm has been already devised. 
Quantum spin models can also be dealt with in the 
present scheme. In a loop-cluster algorithm Jlq , the 
partition function is expressed as a sum of classical (non- 
quantum) weight over spin configurations and graphs. 
The graph degrees of freedom can be divided into a con- 
tinuous part (the locations of the graph elements in the 
imaginary time axis) and a discrete part (the number and 
the types of the graph elements). The present scheme 
can be applicable to the latter discrete part of graph de- 
grees of freedom. The work in this direction is now under 
progress and will be reported elsewhere |ll| . 
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APPENDIX A: RELATIONSHIP BETWEEN TWO 
DENSITIES OF STATES 

We derive the exact relationship between g(E) and 
Cl(n b ). We first define a parameter x as; 

x = exp(l/T). 
Using this x, we can express the partition function as 



Z{x)= 9(E) 



„-E 



X 



E=-N P 

In terms of the number of bonds, it is written as 

JV P 



z( x ) = fi M 0* - !) r 



n b =0 

The range of the energy is — Np < E < and that of the 
number of bonds is < n b < Np where Np is the total 
number of nearest neighbor pairs of spins. Differentiating 
the above two equations I times with respect to x, we 
obtain 



J*-Z(x) = llg(-l) 



1! 



g(-l-l)x 



(1 + 2) 



2! 

-Z(x) = + 



g(-l -2)x 2 + 



([+!)! 
1! 



fi(Z + l)(a;-l) 



(1 + 2)1 
2! 



Q(l + 2) (x- 1) 2 + 



(4) 



By comparing these two equations after taking the limit 
T — > oo, (x — > 1), we arrive at the relation of O (Z ) and 
5(0: 



ni! 

or 

N P -n b 

0(n 6 )= Y, 1 "" ' 

3=0 



g(-n b -j). (0<n b <N P ) 



(5) 



By setting x = in Eq. (0), we obtain 

dN+E , 

9(e) = E (-iy(- E 1 +J 

3=0 V 

(—Np < E < 0). 



fi(-B + j). 



(0) 



Equation (||) is useful for obtaining Ct(nb) from such 
as those obtained by Beale |l4|]. However, computing 
g(E) from fl(rib) using equation (||) is not practical when 
the estimates of fl(ni,) contains statistical error, because 
the (—iy factor in equation (^) magnifies the relative 
magnitude of the errors. 

Using Eq. (||), we can obtain, for example, the expres- 
sion for the ground state entropy, 



Sq 



g(-N P ) = Q(N P ). 



It should be remarked that the direct outcome of the ac- 
tual simulation is not Q(rib) itself but the relative magni- 
tude of Q(n&)'s. Therefore, in order to obtain an estimate 
of Q(Np), we have to use the fact that 



n(o) = q n 

for the Q state Potts model. With this equation, the 
absolute magnitude of f2(rib) can be determined. In other 
words, if fl(Np) is the direct outcome of the simulation 
and therefore proportional to Q(Np), the entropy is given 
by 



,S _ 



n(N P 

f2(0) 



x Q 



N 



APPENDIX B: THE JANKE-KAPPLER 
ALGORITHM 

Here, our implementation of the Janke-Kappler algo- 
rithm 1 12 1 is described. For a given spin configuration 
and a graph, we start with making a random choice of 
a nearest neighbor pair of sites. With some probability, 
we remove the bond if there is one already on the chosen 
pair, whereas we place a new bond if there is no bond on 
the pair and if the pair is satisfied, again probabilistically. 
We say a pair (i,j) is satisfied if Oi = Uj. In cither case, 
the probability for updating is of the heat-bath type; 



P(G'\S,G) ee 



W(S, G') 



W{S, G) + W(S, G') 



where G' is the graph in the proposed final state. 
W(S, G) is defined as 

W(S, G) ee V(G)A(S,G). 

where V(G) is the trial weight which is adoptively ad- 
justed. To be more specific, if there is a bond already on 
the chosen pair, we remove it with probability 

V(n b - 1) 



V(n b ) + V(n b -1)' 



If there is no bond and if the pair is satisfied, we place a 
new bond to the pair with probability 

V(n b + 1) 



V(n b )+V(n b + l)- 
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If the pair is not satisfied, we leave it unconnected. We 
repeat this procedure many times so that every nearest 
neighbor pair is chosen and examined once on the aver- 
age. After these repetitions, we 'flip' all the clusters of 
sites with probability one half. 

If V(rib) is simply written as v nb with some constant v, 
as is the case with the original weight Vq, the decisions of 
placing bonds can be made for each nearest neighbor pair 
independently. In such a case, the resulting algorithm is 
nothing but the Swendsen-Wang algorithm p0|| . How- 
ever, since the adoptively chosen V(n},) is not in general 
factorized, the decisions are dependent. Therefore, in 
the original Janke and Kappler's method only one near- 
est neighbor pair is examined at a time. For this reason 
the graph in the Janke-Kappler algorithm can change 
only gradually. In general, this is a disadvantage because 
there may be some unfavorable region or a 'barrier' in 
the bond-number space which hinders the random walker 
from moving from one side of it to the other. This dis- 
advantage can be removed by the modification proposed 
in the main text, in which the random walker can jump 
from one side to the other in one step without hitting the 
barrier. 

In a practical applications, we perform some number 
of sweeps to obtain a histogram of rib, H(rib). Then, we 
adjust V(n b ) by 

V(n b ) <= V(n b )/H(n b ). 

With this new weight, we redo the simulation. The whole 
procedure is repeated until H{nb) becomes sufficiently n b 
independent. 
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